function ODraw = TreatImage(rawimages,type_img,handles,expp)

% Size of hdf5 image group
[rows,cols,files] = size(rawimages); % rows = horiz. pixels; cols = vert. pixels

%%% TODO: Check validity of image, has pixels, correct size, has been
%%% binned etc.
if rows*cols == 0
    error('No image in HDF file (returned empty matrix)')
end

if  type_img == 3 % Absorption
    if handles.useBackground
        elec = rawimages(:,:,3);
    else
        elec = 0;
    end
    
    if handles.scaleBackground
        sampleIndX = 1:200; % Atoms fall in +x direction so there should never be atoms in this region.
        sampleIndY = 1:200;
        scaleS = mean(mean(rawimages(sampleIndX,sampleIndY,1)));
        scaleB = mean(mean(rawimages(sampleIndX,sampleIndY,2)));
    else
        scaleS = 1;
        scaleB = 1;
    end
    
    num = rawimages(:,:,1) - elec;
    den = rawimages(:,:,2)*scaleS/scaleB - elec;
    
    %scaling the baselevel of the atom image
    numforscale=10;  %using a 10x10 pixel area to determine scaling fraction
    numleftup = sum(sum(num(1:numforscale,1:numforscale)));
    numleftdn = sum(sum(num(rows-numforscale:rows,1:numforscale)));
    numrightup = sum(sum(num(1:numforscale,cols-numforscale:cols)));
    numrightdn = sum(sum(num(rows-numforscale:rows,cols-numforscale:cols)));
    baselevelnum = numleftup + numleftdn + numrightup + numrightdn;
    
    denleftup = sum(sum(den(1:numforscale,1:numforscale)));
    denleftdn = sum(sum(den(rows-numforscale:rows,1:numforscale)));
    denrightup = sum(sum(den(1:numforscale,cols-numforscale:cols)));
    denrightdn = sum(sum(den(rows-numforscale:rows,cols-numforscale:cols)));
    baselevelden = denleftup + denleftdn + denrightup + denrightdn;
    blscale = (baselevelden/baselevelnum); %resulting baselevel scaling
    modnum=num*blscale;
    
    den(den==0) = 1; % Remove zeros from denominator
    
    ODraw = -log((num./den)');  % MATRIX ROTATED HERE
    
    ODrawmod = -log(modnum./den);
    ODraw(isinf(ODraw)) = 0;
    ODraw(isnan(ODraw)) = 0;
    ODraw = real(ODraw);
    
elseif type_img == 1 % Fluorescence
    ODraw = rawimages';         % MATRIX ROTATED HERE
    ODraw = ODraw - min(min(ODraw));
else
    error('Can''t determine imaging method');
end

if handles.doFilterNoise
    % Removes rings from the FFT of the image to remove noise
    Ny = 1:expp.short_px_nb;   %HORIZONTAL DIMENSION OF EXPERIMENT
    Nx = 1:expp.long_px_nb;

    ODf = fftshift(fft2(ODraw));
    
    x = linspace(0,1,length(Nx));
    y = linspace(0,1,length(Ny));
    [X,Y] = meshgrid(x,y);
    
    r1 = handles.innerFFTrad;
    r2 = handles.outerFFTrad;
    c1 = ((X-0.5).^2 + (Y-0.5).^2) > r1.^2;
    c2 = ((X-0.5).^2 + (Y-0.5).^2) < r2.^2;
    mask = find(c1 & c2);
    
    ODf(mask) = 0;
    
    ODraw = real(ifft2(ifftshift(ODf)));
end

if handles.doFilterDC
    % Removes rings from the FFT of the image to remove noise
    Ny = 1:expp.short_px_nb;   %HORIZONTAL DIMENSION OF EXPERIMENT
    Nx = 1:expp.long_px_nb;

    ODf = fftshift(fft2(ODraw));
    
    x = linspace(0,1,length(Nx));
    y = linspace(0,1,length(Ny));
    [X,Y] = meshgrid(x,y);
    
    r = handles.dcCutoff;
    c = ((X-0.5).^2 + (Y-0.5).^2) < r.^2;
    mask = find(c);
    
    ODf(mask) = 0;
    
    ODraw = real(ifft2(ifftshift(ODf)));
end

end